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Abstract. We derive Markovian master equations for single and interacting 
harmonic systems in different scenarios, including strong internal coupling. By 
comparing the dynamics resulting from the corresponding master equations with 
numerical simulations of the global system's evolution, we delimit their validity regimes 
and assess the robustness of the assumptions usually made in the process of deriving 
the reduced Markovian dynamics. The results of these illustrative examples serve to 
clarify the general properties of other open quantum system scenarios subject to be 
treated within a Markovian approximation. 
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1. Introduction 

It is widely assumed that one of the crucial tasks currently facing quantum theorists 
is to understand and characterize the behaviour of reahstic quantum systems. In 
any experiment, a quantum system is subject to noise and decoherence due to the 
unavoidable interaction with its surroundings. The theory of open quantum systems 
aims at developing a general framework to analyze the dynamical behaviour of systems 
that, as a result of their coupling with environmental degrees of freedom, will no 
longer evolve unitarily. If no assumptions are made concerning the strength of the 
system-environment interaction and the time-correlation properties of the environment, 
the dynamical problem may become intractable, despite that the functional forms of 
very general evolutions can be derived [1]. However, there exists a broad range of 
systems of practical interest, mostly in quantum optics and in the solid state physics, 
where it is possible to account for the observed dynamics by means of a differential 
equation for the open system's density matrix derived in the context of Markovian 
processes. Such a differential equation, the so-called Markovian (or Kossakowski- 
Lindblad) master equation, is required to fulfill several consistency properties such as 
being trace preserving and satisfying complete positivity [21 El IH El El El El HO] . 

However, from the theoretical point of view, the conditions under which these type 
of equations are derived are not always entirely clear, as they generally involve informal 
approximations motivated by a variety of microscopic models. This leaves open the range 
of validity of these equations, and which in some circumstances can lead to non physical 
evolutions. The situation becomes even worst as the complexity of the open system 
increases. In particular, it is not an easy question to decide whether the dynamics of a 
composite, possibly driven, quantum system can be described via a Markovian master 
equation, and if so, in what parameter regime. Actually, several groups have recently put 
forward operational criteria to check for deviations from Markovianity of real quantum 
evolutions [IH [121 [131 [H] . 

The main propose of this work is to study such interacting open quantum systems, 
and show that there are Markovian master equations close to the real dynamics, 
characterizing the range of validity of each one. To this aim we have chosen a system 
consisting of quantum harmonic oscillators, as one can easily follow the exact dynamics 
using numerical simulations of a particular, but wide class of simple states, the so-called 
Gaussian states. Moreover, the proposed method is general enough to be applicable 
to non-harmonic systems and, in particular, when the coupling between oscillators is 
sufficiently weak so that their local dynamics is effectively two-dimensional, we expect 
the conditions obtained for strict Markovianity to be directly applicable to systems of 
interacting qubits. 

The damped harmonic oscillator is the canonical example used in most references 
to discuss both Markovian and non-Markovian open system dynamics (see for instance 
[HI [m [ini [171 [El [El [20] and references therein) and exact solutions in the presence of 
a general environment are known [21j. The dynamics of coupled damped oscillators. 
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including those interacting with a semiclassical field, are significantly less studied, 
with most analysis focusing on evaluating the decoherence of initially entangled states 
provided that certain dynamical evolution, Markovian or not, is valid |22]. Recently, an 
exact master equation for two interacting harmonic oscillators subject to a global general 
environment was derived |23] . Here we will focus on the derivation of Markovian master 
equations for interacting systems. We will focus on a scenario where two harmonic 
systems are subject to independent reservoirs and present a detailed study based on the 
numerical simulation of the exact dynamics. The advantage of this approach is that 
it allows us to compute not only quantities for the damped system but also for the 
environment. This enables us to check the rigour of some of the assumptions usually 
made in obtaining a Markovian master equation and assess their domain of validity. 

We have extensively studied three damped systems. For completeness, we start 
our analysis by considering a single harmonic oscillator (section [2]) and subsequently 
move to the core of our study by analyzing the dynamics of two interacting harmonic 
oscillators (section [3]), finding Markovian master equations for both weak and strong 
internal coupling. We finally address the dynamics of an harmonic oscillator driven 
by a semiclassical field (section H]), where different Markovian master equations have 
been obtained and studied depending on the values of the external Rabi frequency and 
the detuning from the oscillator's natural frequency. To make the reading more fluent, 
details of the simulations and the derivation procedure are left for the appendices. 

In the following two introductory sub-sections, and with the aim of setting up the 
notation and making the presentation as self contained as possible, we present a brief 
discussion of how Markovian master equations are obtained in the weak coupling limit 
(section II. ip . and present a short review of the properties of the harmonic oscillator 
Gaussian states, which will be used in subsequent sections (section [L2|) . 

1.1. Markovian Master Equations 

To derive Markovian master equations we follow the approach of projection operators 
initiated by Nakajima [24j and Zwanzig [25], see also [Sj [151 [E] for instance. In 
this method we deflne in the Hilbert space of the combined system and environment 
= Tis <^T-Le two orthogonal projection operators, Vp = Ti e{p) ® Pe and Q = \ — V . 
Here p G 23(7/) is the combined state and pe G '^{1-Le) is a flxed state of the 
environment, which we choose to be the real initial (thermal, = 1) state, 

Pe = Ah = exp{-HE/T){TT[exp{-HE/T)]}-\ 

Note that Vp gives all the necessary information about the reduced system state ps, so 
to know the dynamics of Vp implies that one knows the time evolution of the reduced 
system. 

We then assume that the dynamics of the whole system is given by the Hamiltonian 
H = Hs+He+oV 1 where Hs and He are the individual Hamiltonians of the system and 
environment respectively and V describes the interaction between them with coupling 
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strength a. Working in the interaction picture {h = 1), 

p{t) = exp[t{Hs + HE)t]p{t) exp[~t{Hs + HE)t], 
and analogously for V(t), we obtain the evolution equation 

|p(t) = -z«[\/(t),p(t)]=aV(t)p(t). (1) 

For the class of interactions that we are interested in TTE[V{t)pE] = Tr£;[V^(t)pth] = 
0, which implies 

VV{t)V = 0, (2) 

as can be easily checked by applying it over an arbitrary state p G 23(7/). It is not 
difficult to redefine the interaction Hamiltonian such that this always holds, see for 
example [T0l|20]. 

Our aim is to obtain a time-evolution equation for Vp under some approximation, 
in such a way that it describes a quantum Markovian process. To this end, we apply 
the projection operators on equation ([1]), introducing the identity 1 = V + Q between 
V{t) and p(t), 

= aVV{t)Vp{t) + aVV{t)Qp{t), (3) 

^Qp(t) = aQV(t)Vp(t) + aQV{t)Qp{t). (4) 
The solution of the second equation can be written formally as 

Qpit) = g{t,to)Qp(to) + a [ dsg{t,s)QVis)Vp{s). 

J to 

This is nothing but the operational version of the variation of parameters formula for 
ordinary differential equations (see for example [26], [27]), where the solution to the 
homogeneous equation 

|Qp(t) = «QV(t)Qp(t) 

is given by the propagator 

^(t,,) = re"/>*'W), 

where T is the time-ordering operator. Inserting the formal solution for Qp{t) in 
yields 

Ipp(t) = aVV{t)Vp{t) + aVV{t)g{t,to)Qp{to) 
at 

W [ dsVV{t)g{t,s)QV{s)Pp{s). 

J to 

We now assume that the initial state of the system and bath are uncorrelated, so that 
the total density operator is factorised into p(to) = ps(^o) ® Pth- From this we find 
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Qpito) = 0, which was guaranteed by our choice of V as projecting onto the initial 
state, and then by using (2) we finally arrive at 

^Vp{t)= f dslCit,s)Vp{s), (5) 



dt 



to 



with kernel 



)C{t,s) = a^VV{t)g{t,s)QV{s)V. 



Equation ([5]) is still exact. We now consider the weak coupling limit, by taking the 
kernel at lowest order in a, 

/C(t, s) = a^VV{t)QVis)V + C(a^), (6) 

so that by again using condition ([2]) we get a Born approximation for (jS]): 



dt 

which implies 



d /■* 

'^Vp{t) = a^ dsVV{t)V{s)Vp{s), 



to 



j^psit) = -a^ dsTiE[V{t), [V{s)rps{s)®p,^]]. (7) 

Note that we are not asserting here that the state of the bath is always pth? the term 
Ps{s) ® pth appears just as a result of the application of the projection operator (see 
discussion in section I2.3.3p . Now we take the initial time to = and an elementary 
change of variable s by t — s in the integral yields 

-ps{t) = -a" j ds TiEiVit), [V{t - s),ps{t -s)(g) pth]]. 

We expect this equation to be valid in the limit a — )■ 0, but in such a limit the change 
in ps becomes smaller and smaller and so if we want to see dynamics we need to rescale 
the time by a factor [21 SJ |5] otherwise the right side of the above equation goes 
to zero. Thus in the limit a — )■ the integration is extended to infinity. However in 
order to get a finite value for the integral, the functions TYE[V(t), [V(t — s),pb]] must 
decrease appropriately. In particular this implies that they should not be periodic, which 
requires that the number of degrees of freedom in the environment must be infinite, as 
otherwise there will be a finite recurrence time. Moreover, as ps changes very slowly 
in the limit a — )■ 0, we can take it as a constant inside width tb around s = where 
Tr^[K(t), [V(t — s),pb]] is not zero, and so finally we obtain 
d ^ ^ 

-ps{t) = -a'J^ dsTTE[V{t),[V{t-s),ps{t)®p,^]]. (8) 

These informal arguments contain the basic ideas behind the rigorous results obtained 
by Davies [HIS]. 

Since we have started from a product state p(to) = Ps(^o) ® Pth, we require, for 
consistency, that our evolution equation generates completely positive dynamics. The 
last equation does not yet warrant complete positivity in the evolution [8] , and so we need 
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to perform one final approximation. To this end, note that the interaction Hamiltonian 
may be written as: 

V = Y,Ak®Bk, (9) 

k 

where each can be decomposed as a sum of eigenoperators of the superoperator 
[Hs. ■] 

Ak = Y,M^). (10) 

where 

[Hs.Ak{p)] = -uAk{p). (11) 

This kind of decomposition can always be made [3l [10]. On the other hand, by taking 
the Hermitian conjugate, 

[Hs.A\{u)] = uAl{ul 

and since V is self-adjoint, in the interaction picture one has 

V{t) = Y^e-^^'A,{v)®B,{t) = $^e-*4(z.) ®5t(t). 

Now, substituting the decomposition in terms of Afc(z/) for V(t — s) and AI^u) for V(t) 
into equation ([8]) gives, after expanding the double commutator, 

j/s{t) = $^$^e^('''-'^)*r,,(z.)[A,(z/)p5(t),4(^')] 

uy k,i 

+e^(^-'^')*ri,(z.)[A,(z/'),P5(t)4(^)], (12) 
where we have introduced the quantities 

TkM = a" / cise^'^^Tr Bl{t)B,{t - s)p,^ 
Jo L 



Blis)B,ptt\ , (13) 

with the last step being justified because pth commutes with exp{iHEt). 

In equation ( !T2|) the terms with different frequencies will oscillate rapidly around 
zero as long as |z/' — z/| ^ a^, so in the weak coupling limit these terms vanish to obtain 

j/sit) = $^$^r,X^)[A,(^)p5(t),4(^)] +n,.(^)[A,(z.),p^(t)4(z.)]. (14) 

u k,l 

Now we decompose the matrices Tk,i{v) as a of sum Hermitian and anti-Hermitian parts 
where the coefficients 
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and 



ki 



form Hermitian matrices. In terms of these quantities (031) becomes 



j^Psit) = -i[Hi^srPs{t)] + V[Ps{t)], 



where 



u kl 

is a Hermitian operator which commutes with Hs-, as a consequence of (fTTj) . This is 
usually called the Shift Hamiltonian, since it produces a renormalization of the free 
energy levels of the system induced by the interaction with the environment. The 
dissipator is given by 



it, 



A,{v)ps{t)Al{v) - ]^{Al{v)A,{v)rPsit)} 



V k,e 

Returning to Schrodinger picture, the time-evolution equation is then just 
d 



dt 



psit) = -i[Hs + i/LS, Ps{t)] + V[ps{t)]. 



(15) 



Note that the matrices jk/i^) are positive semidefinite for every u, this is a consequence 
of the Bochner's theorem [28] , that is, it is easy to check that the correlation functions 



Tr 



Bl{s)B,p,^ 



of them. Wit 



are functions of positive type, and '^k/{i^) are just the Fourier transform 
1 this final remark we conclude that the equation ( |T5l) generates a 
completely positive semigroup [6] and so defines a proper Markovian master equation, 
i.e. a completely positive semigroup. 



1.2. Gaussian States 

We saw in the last section that to avoid a finite recurrence time, the number of 
environment degrees of freedom should strictly tend to infinity. However, in practice, the 
recurrence time grows very rapidly with the size of the environment and so one can still 
test the validity of such equations with only a finite, yet still large environment model, 
as long as the domain of interest is restricted to early times. The prototypical example 
of which is afforded by a collection on n harmonic oscillators. In fact, such models are 
often explicitly included in master equation derivations both due to their easy handling 
and due to realistic physical justification. Phenomenologically speaking, they correctly 
describe both quantum Brownian motion and the derivation of Langevin style equations 
from first principles [16]. However, they also provide a convenient numerical testing 
ground as the number of variables needed to model such systems scales polynomially in 
the number of degrees of freedom. This is because the harmonic oscillator falls into a 
class of quantum states known as Gaussian states, which are entirely characterised by 
their first and second moments. We now review some of their basic properties 
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For any system of n canonical degrees of freedom, such as n harmonic oscillators, or 
n modes of a field, we can combine the 2n conjugate operators corresponding to position 
and momentum into a convenient row vector, 



R= {Xi,X2,...,Xn,PuP2,---,PnV- 

The usual canonical commutation relations (CCR) then take the form 

[Rk, Rl] = iho'kl, 

where the skew-symmetric real 2n x 2n matrix a is called the symplectic matrix, 
the choice of R above, a is given by. 



a 





-1, 



In 





One may also choose a mode-wise ordering of the operators, R = (xi,pi, ...,Xn,i^n 
which case the symplectic matrix takes on the form. 



(16) 

(17) 
For 

(18) 



a 



e 







1 





(19) 



Canonical transformations of the vectors S : R ^ R' are then the real 2n— dimensional 
matrices 5* which preserve the kinematic relations specified by the CCR. That is, the 
elements transform as R' = SabRb, under the restriction. 



SaS' 



a. 



(20) 



This condition defines the real 2n-dimensional symplectic group Sp(2n, R). For any 
element S G Sp(2n, R), the transformations —S, and are also symplectic 
matrices, and the inverse can be found from 5*"^ = aS"^a~^. The phase space then 
adopts the structure of a symplectic vector space, where f|T9l) expresses the associated 
symplectic form. Rather than considering unitary operators acting on density matrices 
in a Hilbert space, we can instead think of all the quantum dynamics taking place on 
the symplectic vector space. Quantum states are then represented by functions defined 
on phase space, the choice of which is not unique, and common examples include the 
Wigner function, Q-function and the P-function. Often one has a particular benefit 
for a given physical problem, however for our purposes we shall consider the (Wigner) 
characteristic function XpiO^ which we define through the Weyl operator 



^ e R"' 



2n 



(21) 



as 



X,iO = npW^]- (22) 

Each characteristic function uniquely determines a quantum state. These are related 
through a Fourier- Weyl transform, and so the state p can be obtained as 

^ '' d'-^Xpi-OW^. (23) 



P 



{2< 



2n 
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We then define the set of Gaussian states as those with Gaussian characteristic functions. 
Equivalent definitions based on other phase space functions also exist, but for our choice 
we consider characteristic functions of the form, 

= Xp(0)e-^«"^«+^"^ (24) 

where C is a 2n x 2n real matrix and D G R^" is a vector. Thus, a Gaussian characteristic 
function, and therefore any Gaussian state, can be completely specified by 2n^ + 3n real 
parameters. The first moments give the expectation values of the canonical coordinates 
dj = TT[Rjp\ and are related to D by d = a~^D, while the second moments make up 
the covariance matrix defined by 

Cj,k = 2Re TMRj - {Rj)p){Rk - {Rk)p)]. (25) 

These are related to C by the relation € = a^Ccr. It is often the case that only the 
entanglement properties of a given state are of interest. As the vector d can be made 
zero by local translations in phase space, one can specify the state entirely using the 
simpler relation, 

Cj-fc = 2Re Tr[pi?ji?fc]. (26) 

However, in this work we shall predominantly use the relation (1251) . Using this 
convention, we mention two states of particular interest; the vacuum state, and the 
n-mode thermal state. Both take on a convenient diagonal form. In case of the vacuum 
this is simply the identity C = l2n, while for the thermal state the elements are given 
by 

where Uj is the frequency of the j^^ mode, and the equilibrium temperature is given by 
T. 



1.2.1. Operations on Gaussian States We now consider Gaussian transformations. 
As the Rj are hermitian and irreducible, given any real symplectic transform S, the 
Stone- Von Neumann theorem tells us there exists a unique unitary transformation Us 
acting on Ti such that UsW^Ul = Ws^. Of particular interest are those operators, 
Ug, which transform Gaussian states to Gaussian states. To this end, we consider the 
infinitesimal generators G, of Gaussian unitaries Ug = e~"*^ = 1 — ieG + C(e^). Then 
to preserve the (Weyl) canonical commutation relations, the generators G must have 
the form G = Yl^l:=i9jk{RjRk — RkRj)/'^ I2j- It follows that Hamiltonians quadratic 
in the canonical position and momentum operators (and correspondingly the creation 
and annihilation operators) will be Gaussian preserving, in particular, the Hamiltonian 
for n simple harmonic oscillators, H = Ylj=i^j^]^j- It is for this reason that harmonic 
oscillators provide such a useful testing ground for many body systems. 

An additional, though simple, property worth highlighting is the action of the 
partial trace. Using the expression for the density matrix (1231) . it is straightforward to 
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see the effect of the partial trace operation on the characteristic function. If we take a 
mode-wise ordering of the vector R = {Ri,R2), where Ri and R2 spht two subspaces 
of ni and ^2 conjugate variables corresponding to partitions of the state space of p into 
"H = "Hi ® 'H2, then the partial trace over is given by 



f28) 



That is, we need only consider the characteristic function associated to the vector 
Ri. At the level of covariance matrices, we simply discard elements corresponding to 
variances including any operators in R2, and so the partial trace of a Gaussian state 
will itself remain Gaussian. 

Finally, we make some remarks regarding closeness of two Gaussian states. Given 
pi and p2 the fidelity between them is defined as F{pi, P2) = (Tr a/ , and is a 
measure of how close both quantum system are each other. Actually a distance measure 
can be defined as Db = y/l — F which is essentially the same as the Bures distance [29] 
(^-^Bures(Pi' P2) = 2 — 2 a/-F(pi , P2) j • This distance will be very useful for quantifying 
how well the dynamics generated by a Markovian master equation approximate the real 
one. 

In general the fidelity is quite difficult to compute, however in the case of Gaussian 
states Scutaru has given closed formulas in terms of the covariance matrix |30]. For 
example, in case of one mode Gaussian states pci and pg2, with covariance matrices 
C^^^ and C-^-* and displacement vectors and d^"^^ respectively, their fidelity is given 
by the formula 

F(pai,PG2) = ^"^P \-^^ {C^'^ + C^'^y' 5] , (29) 

V A + $ - V $ ^ ^ 

where A = det [C^^^ + C^^)] , $ = det (C^^) - l) det (C^^) - l) and 6 = (rf^^) - rf^^)) . 



2. Damped Harmonic Oscillator 

We will first consider a single harmonic oscillator damped by an environment consisting 
of M oscillators (see figured]). We want to know under which conditions the Markovian 
master equation that we derived in the previous section for the evolution of the damped 
oscillator is valid. To this aim we will approach the exact dynamical equations of the 
whole system when M is large; these will be solved via computer simulation, and we 
can then compare this solution with the one obtained using a master equation. 
The Hamiltonian for the whole system will be given by {h = 1) 

M M 

Wj-ajaj + 5fj(a^aj + aaj). (30) 
i=i i=i 

Note that the coupling to the bath has been considered in the rotating wave 
approximation (RWA), which is a good description of the real dynamics for small 
damping Q ^ max{gj,j = 1, . . . , M} (e.g. in the weak coupling limit) |3T]. 
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Figure 1. Model for a damped harmonic oscillator. The central grey sphere represents 
the damped oscillator which is coupled to a large number of environmental oscillators 
(blue spheres) with different frequencies ujj via the coupling constants gj, These are 
chosen in agreement with an Ohmic spectral density pip . 



For definiteness, in this paper we have chosen to distribute the environmental 
oscillators according to an Ohmic spectral density with exponential cut-off. In the 
continuous limit, this has the form [19] 

M 



where a is a constant which modifies the strength of the interaction and oOc is the 
so-called cutoff frequency. Clearly J{u) increases linearly for small values of u, decays 
exponentially for large ones, and has its maximum at a; = Uc- Of course any other choice 
of spectral density could have been taken, but this in turn would require a re-analysis 
of the master equations' range of validity. 

2.1. Exact Solution 

The exact solution of this system can be given in terms of the time-evolution of the 
collection {a,aj} in the Heisenberg picture [17]. From (130|1 we have 

AI 




(31) 




(32) 



icij = [oj, H] = oOjOj + gja, 

and so by writing A = (a, oi, 02, . . . , om) 
expressed as 




iA = WA, 



(34) 
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where W is the matrix 

/ n gi 92 ■■■ Qm \ 

91 i^i 

92 "^2 



W 



(35) 



\ 9m I 
and the solution of the system will be given by 

A(i)=rA(0), r = e-^^*. (36) 

Analogously, the evolution of the creation operator will be 

- iit = WA^ A\t) = TUt(O), = e'^*, (37) 

where = {a^ , a\, . . . , a\_^)^. 

We can also compute the evolution of position and momentum operators X — 
\{A + A^) andP= ^(A-At), 

x{t) = -[TA(o) + rUt(o)] 

2 

= i{T[X(0)+^P(0)]+Tt[X(0)-zP(0)]} 

^TrX{0)-TiP{0), (38) 

and similarly 

P(t)^TiX(0)+TRP(0), (39) 
in these expressions, Tji and Tj are the self- adjoint matrices defined by 

r=T« + >T,^P«%'f = ™^('r'> . (40) 
So, the time-evolution of the vector R — {x,Xi, . . . , xm,P,Pi, ■ ■ ■ ,Pm)'^ will be given by 
P(t) = A^P(O) = || "r^^^m, (41) 

note that the size of M is 2(M + 1) x 2(M + 1). 

Due to the linearity in the couplings in H, an initial (global) Gaussian state po will 
remain Gaussian at all times t, and so we can restrict our attention to the evolution of 
its covariance matrix 

dj = {RiRj + RjRi) - 2{Ri) (Rj) . (42) 

Particularly, since we are interested in just the first oscillator, we only need the evolution 
of the 2x2 submatrix {Cij;i,j — 1,M + 2}. The evolution of pairs of position and 
momentum operators is 

{Ri{t)Rj{t)) = J2■M^,kMJ,e{M0)M0)): (43) 
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and similarly for products of expectation values {Ri(t)){Rj{t)). So the elements of the 
covariance matrix at time t will be 

C,j(t) = {Ri{t)Rj{t) + Rj{t)Ri{t)) - 2{Ri{t)){R,{t)) 

= Y^Mi,kMjA{Rk{0)Re{0) + ReiO)Rk{0))-2{Rkm{Ri{0))] 

= ^Mi,kMj,iCk,i{0), 

k,e 

and for the first oscillator we have 

Ci,i(t) = J2Mi,kMi,iCkA0) = {Mi,CMi), (44) 

k,e 

Cl,M+2(t) = CM+2,l{t) = ^ All,fcA^A/+2A,f(0) = {Mi,CMm+2), (45) 

k,e 

(t) = MM+2,kMM+2,M0) = {Mm+2, CMm+2), (46) 

k/ 

here (■, ■) denotes the scalar product, and the vectors Aii and A^a/+2 are given by 

Ml = (7Wi,i, A^i,2, . . . , -Mi,2M+2)^, (47) 

■M.M+2 = (A^M+2,1, A^Af+2,2, • • • , A/+2,2M+2)'^- (48) 

More details of how this exact solution is simulated in order to approach the 



Markovian master equation description are given in Appendix A 



2.2. Markovian Master Equation 

The damped harmonic oscillator is a standard example for the derivation of master 
equations (see for example [31 [171 UHl En] ) . The Markovian master equation (flSil is given 

by 

^p(t) = -in[a^a, p{t)] + 7(n + 1) {2ap{t)a^ - a^ap{t) - p{t)a^a) 

+7^ [2a^p{t)a - aa^ p{t) - p{t)aa^) , (49) 
where is a renormalized oscillator energy arising for the coupling to the environment 

n = ^] + A, A = P.V. / rfw-^^, (50) 

(here P.V. denotes the Cauchy principal value of the integral), n is the mean number 
of bath quanta with frequency f2, given by the Bose-Einstein distribution 

-1 



n = nsi^, T) 



exp I p ) - 1 



(51) 



and 7 is the decay rate, which is related to the spectral density of the bath J{uj) = 
T.j9p{^j-^) via 

7 = TiJin). (52) 
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Note that the shift A is independent of the temperature, and although its effect is typical 
small (e.g. [3l [18]) we will not neglect it in our study. For an ohmic spectral density the 
frequency shift is 



In addition, note that the equation (149|) is Gaussian preseving [32], as it is the limit 
of a linear interaction with an environment and so the total system remains Gaussian 
while the partial trace also preserves Gaussianity. 

2. 3. Study of the Approximations 

As a first step, we have plotted the variance of the x coordinate for two different initial 
states of the system, these are a thermal and a squeezed state, see figure |2j The last 
plot clearly illustrates the closeness of the results for the Markovian master equation, 
when compared to the effect of the Lamb shift. To explore this further, we now study 
several effects which pertain to the validity of this equation, by calculating the distance 
(in terms of the fidelity) between the simulated state p^g' and the state generated by 
the Markovian master equation p^"*. 

2.3.1. Discreteness of the hath Due to the finite number of oscillators in the bath, we 
can only simulate inside a bounded time scale free of the back-action of the bath. This 
produces revivals in the visualized dynamical quantities for times t < tr, where tr is the 
recurrence time of the bath. Of course, the time after which these revivals arise increases 
with the number of oscillators in the bath, and roughly speaking it scales as tr oc M. 
This behaviour is shown in figure [3l where the distance between the simulation and the 
Markovian master equation for a system initially in a thermal state with temperature 
Ts = 30 is plotted as a function of the time and the number of oscillators. 

2.3.2. Temperature It is sometimes claimed that for ohmic spectral densities the 
Markovian master equation ( l49l) is not valid at low temperatures [TH [19]. Of course, 
one must make clear the context in which this claim is made, and so for definiteness, 
let us focus on the validity with respect to the bath temperature. A detailed discussion 
of this situation can be found in the book by Carmichael [18]. There the argument 
is based on the width of the correlation function Ci2(t) = TT[Bi{s)B2Pth], where 

-_igjaj, which increases for an Ohmic spectral density as the bath 
temperature decreases. More specifically, in the derivation of the Markovian master 
equation two kinds of correlation functions appear. 




where Ei is the exponential integral function defined as 




Ci2(s) = Tr[5i(s)52Pth] = 5^^/fc^/,e*"^-Tr[a]afePth] 
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System in a Themal State Ts = 30 
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System in a Squeezed State (Aj;)^ = [2(A;))]"^ = 5/2 
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Figure 2. Comparison of the evolution of 2(Aa;)^ for an initially thermal and squeezed 
(vacuum) state. The bottom plot shows the effect of the Lamb shift, which produce a 
"slippage" in the squeezed state variances. 




Time 



Figure 3. Color map showing the dependency of the recurrence times with the size 
of the bath. The rest of the parameters are the same as in figured 
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M 
3 

and 

C2i{s) = Tr[52(s)5ipth] = 5^(?,(?,e-*-^-Tr[a,4pth] 

M 

j 

We may call Ci2(s) = C(— s,T) and 6*21(5) = C{s,T) + Co(s), and so in the continuous 
limit 

POO POO ,2 

Co{s) = / J{u)e-''^'du = a ue-''^^'-'^^'^du 



[isujc + 1 



>2' 



and 



T) = J{uj)e-'^'n{uj, T)duj = aT\ (^2, 1 - isT + , 

where here ({z, q) = X^fcLo [(q+k)2]=^/2 so-called Hurwitz Zeta function, which is a 

generalization of the Riemann zeta function ({z) = ({z, 1) j33] . 




s Tb 

Figure 4. On the left, the absolute value of the correlation function is plotted for 
several temperatures while the FWHH as a function of temperature is represented on 
the right. 

In the left plot of figure HI the absolute value of C{s,T) is plotted for different 
temperatures. Note that the spreading of the correlation function is mainly caused 
by its "height" decrease, that is, in the limit T — )■ 0, C(s,T) — 0. So one may also 
expect that the contribution of these correlations to the motion becomes less important 
as T — )■ 0, in such a way that the problem of the infinite width can be counteracted. 
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and this is indeed what seems to happen. To visuahze this more carefully we have 
plotted in the right of figure HJthe full weight at half height (FWHH) for both Cq{s) and 
C{s,T). In order to make valid the Markovian approximation, the typical time scale 
for the evolution of the system due to its interaction with the bath ts must be large in 
comparison with the decay time tb of the correlation functions. Loosely speaking, this 
can be characterized by the FWHH. 

From figure Hlone sees that for small temperatures tb (i.e. FWHH) is quite large, so 
it is expected that the Markovian approximation breaks down for values of T such that 
Ts ^ Tb- However if a is small enough this will happen for values where the contribution 
of C{s,T) to the convolution integrals is negligible in comparison with the contribution 
of Cq{s), whose FWHH will remain constant and small with respect to ts- As a rough 
estimation, using the parameters in figure [2], we find that to get a value of the FWHH 
comparable with ts ~ l/\/a ~ 22.4, we need a temperature of at least T ~ 0.05. Both 
contributions enter in the Markovian master equation derivation via some convolution 
with the quantum state and one oscillating factor. We may get a very informal idea 
of how both contributions matter by looking at their maximum values at s = 0, for 
example C{s = 0,T = 0.05) = 3.27391 x 10"^ and Co{s = 0) = 0.018, and so it is clear 
that C{s,T = 0.05) will not have a large effect on the dynamics. For large temperatures 
the FWHH of C{s,T) remains small though now larger than Cq{s), so it is expected 
that in the limit of high temperatures the accuracy of the Markovian master equation 
stabilizes to a value only a little worse than for T = 0. 

All of these conclusions are illustrated in figure [5], where the fidelity between the 
state from the simulation and that from the Markovian master equation is plotted. The 
behaviour at very early times is mainly related to the choice of the initial state of the 
system, and reflects how it adjusts to the state of the bath under the Markovian evolution 
[3^ . different tendencies have been founded depending on the choice of initial state. 
However the behaviour with temperature is visible at longer times (since tb ~ FHWW 
increases with T) which is in agreement with the conclusions drawn from the correlation 
functions (see small subplot). At zero temperature (blue line) the results are in closest 
agreement, however, as the temperature is increased to T = 0.1 the correlation function 
broadens, which leads to a degradation (albeit small) in the modelling precision. As the 
temperature increases further, the influence of this correlation function becomes more 
important and the FWHH decreases to a limiting value (see the plot on the right of 
figure Hj), this convergence is reflected by the red, cyan and purple lines which show that 
the accuracy at large temperatures stabilizes to only a little worse than that at T = 0, 
as was expected from figure HI 

In summary, the Markovian master equation (H9l) does not properly describe the 
stimulated emission/absorption processes (the ones which depend on C{s,T)) for low 
temperatures, however the temperatures when this discrepancy is apparent are so 
small that the contribution from stimulated process are negligible in comparison with 
spontaneous emission, and so the discrepancy with the Markovian master equation is 
never large. 
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Parameters: 
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Figure 5. Fidelity between the simulated state ps and that given by the Markovian 
master equation time evolution pA/ , for several temperatures. For large times (see inset 
plot) temperature does not play a very significant role in the accuracy while for small 
times the accuracy depends mainly on the choice of the initial state of the system (see 
discussion in the text). 



2.3.3. Assumption of factonzed dynamics p{t) = ps(t) (g) pth In the derivation of 
the Markovian master equation, one can arrive at equation ([7]) by iterating the Von- 
Neumann equation ([T]) twice and assuming that the whole state factorizes as p(t) ~ 
Psit) ® Pth at any time ([U [TOl [HI ISO])- This assumption has to be understood as 
an effective model for arriving at equation ([7]) without the use of projection operator 
techniques, however it does not make sense to assume that the physical state of the 
system is really a factorization for all time. Taking advantage of the ability to simulate 
the entire system we have plotted the distance between the simulated whole state p(t) 
and the ansatz psit)<S)pth as a function of time, see figure [61 On the left we have plotted 
the distance for M = 350 oscillators in the bath, actually we have checked from several 
simulations that the results turn out to be independent of the number of oscillators as 
long as the maximum time is less than the recurrence time of the system. From figure 
[3] we see that t = 50 is less than the recurrence time for M = 175, and so we have used 
this value and plotted the distance for different coupling strengths on the right. It is 
clear that this distance is monotonically increasing in time (strictly, in the limit of an 
environment with infinite degrees of freedom), and the slope decreases with coupling 
strength. In section 11.11 we pointed out that the weak coupling approach make sense if 
the coupling is small and the environment has infinite degrees of freedom. This fits with 
the usual argument to take p ~ Ps{t) 8) Pth in niore informal derivation of Markovian 
master equations, that is "the state of the environment is not so affected by the system" , 
but we stress again that this is an effective approach, without any physical meaning on 
the real state p. 
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Figure 6. Distance between tlie simulated states ps{t) (8 pth and p{t) as a function of 
time; on the left side, different numbers of oscillators in the bath are plotted (obtaining 
the same result on that time scale), and on the right side, different values of the coupling 
constant have been taken. 



3. Two Coupled Damped Harmonic Oscillators 

We now consider two coupled harmonic oscillators, which for simplicity we take to have 
the same frequency Qi = Q2 = ^, and each locally damped by their own reservoir (see 
figure [7j), the Hamiltonian of the whole system is 

H = Hqi + Hq2 + Vi2 + Hbi + Hb2 + ViBl + V2B2-, (53) 

where the free Hamiltonians are given by 

Hqi = Qa\ai, Hq2 = ^ala2, 

M M 

Uijaljaij, Hb2 -'^uJ2jalja2j, 
i=i i=i 
with the couplings to the baths, 

M 

ViBi = '^gij{a\aij + aia\j), 
i=i 

M 

V2B2 = ^fi'2j(4o2j + a2alj), 

and the coupling between oscillators, 
V12 = f3{a\a2 + ai4). 

Again we have employed the rotating wave approximation, and so we assume Q ^ p. 
For the case of f2 ~ /3 we must keep the antirotating terms 0102 and ajog. However 
note that the eigenfrequencies of the normal modes become imaginary if u < 2/3 (see 
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for example PJ^J) and the system then becomes unstable, so even when keeping the 
antirotating terms, we must limit /3 if we wish to keep the oscillatory behaviour. 
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Figure 7. The same model as figured] for the case of two damped harmonic oscillators 
coupled together with strength f3. 



3.1. Exact Solution 

For the exact solution, the extension to two oscillators follows closely that of a single 
damped harmonic oscillator. Again, we work in the Heisenberg picture, and wish to solve 
for the vector A= (ai, an, . . . , aiM, 02, 021, • • • , 0'2m)^ ^ given the differential equation, 

iA = WA, (54) 

where W is now given by the matrix 

/ fii gu ■■■ giM (3 \ 

9n UJ12 



W 



9lM 



^2 ^721 
^721 ^21 



92M 



\ g2M ^2IVI J 

The simulation process is then analogous to that of section 12.11 



(55) 



3.2. Markovian Master Equations 

Unfortunately, the derivation of a Markovian master equation for coupled systems 
introduces a number of additional complications. If the oscillators are uncoupled /3 = 0, 
it is obvious that the Markovian master equation for their joint density matrix will be 
a sum of expressions like ( H^ . 

^osit) = -ipa\ai + ^a\a2,ps{t)]+V^[ps{t)]+V2[ps{t)], (56) 
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where 



'jjirij + 1) {2ajps{t)a] - a]ajPs{t) - ps{t)a]a 
+'yjnj (2a^jPsit)aj — aja^jp{t) — ps{t)aja^j 



(57) 



here each frequency shift, decay rate and number of quanta are individually computed 
via equations f l5T]) . f l52|) and fl50l) for each bath j. However for finite intercoupling we 
split the analysis in two subsections. 



3.2.1. Small intercoupling /3 If /3 is sufficiently small to not affect the shift and decay 
rates, one can expect a Markovian master equation of the form 
d 



dt 



psit) = ~t[nalai + nala2 + Vi2,Psit)]+Vi[psit)] + V2[psit)], (5J 



an example of which for coupled subsystems can be found in [36j, and we have given the 



details of a derivation based on projection operators in Appendix B.l In addition, this 
kind of approximation is often made in other contexts such as with damped systems 
driven by a classical field [18] . Such a case will be analyzed in detail in section HI 



3.2.2. Large intercoupling /3 To go further we must work in the interaction picture 
generated by the Hamiltonian Hq = //free + Vu and apply the procedure described in 



section ILTI The details of the derivation are left for Appendix B.2 what is important 
however, is that the non-secular terms oscillate with a phase e^^*'^* so in order to neglect 
them we must impose /3 ^ a, therefore the resultant equation is, in some sense, 
complementary to ( l58ll valid if a > /3. The final Markovian master equation in this 
regime takes the form 



d_ 
di 



Ps(t) = -i[^a\ai + Qala2 + (3 (ai^ + a|a2 ) ,Ps(t)] 



here 



(E) 
k 



2 



jk 



ajps{t)al + ^{aloj, ps{t)} 



a]ps{t)ak + -{ako], ps{t)} 



(59) 



n = n + [Ai{n+) + A2(^^+) + Ai{n_) + A2(^^-)]/4, 

(3 =(3 + [Ai(fi+) + A2(f]+) - Ai(fi_) - A2(fi-)]/4, 
and ii^jf ^ and i^jf are two positive semidefinite Hermitian matrices with coefficients 



K 



(E) 
11 



K 
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K 



(E) 
12 



K 



(E)* 
21 



= {^,in+)[n^in+) + 1] + 72(^]+)[r^2(^^+) + 1] 
+7i(fi_)[ni(fi_) + 1] + ^2{^^)[n2{^-) + l]}/2, 
= {7l(fi+)[nl(^]+) + 1] + 72(^]+)[n2(fi+) + 1] 

-7i(r]_)[r2i(fi_) + 1] - 72(^]-)[n2(^^-) + 1]}/2, 



(60) 
(61) 
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K[f = K^f = [7l(^]+)nl(fi+) +72(fi+)n2(^^+) 

+7i(^^-)ni(fi_) + 72(^^-)^2(^^-)]/2, (62) 
K[f = Kif* = [7i(l^+)ni(l]+) +72(f^+)^2(^^+) 

-7l(^^-)r^l(^^-) - 72(^^-)ri2(^^-)]/2, (63) 

where •jj, Aj and uj are evaluated according to the spectral density and temperature of 
the bath j and Q± = Q ± /3. 

3. 3. Study of the Approximations 

By virtue of the derivation, equations (!58l) and (!59l) preserve both complete positivity 
and Gaussianity (because they arise from a linear interaction with the environment). 
Thus we can test their regimes of validity using simulations of Gaussian states, and 
the appropriate fidelity formulas. In figure M we have plotted the fidelity between both 
states for the Markovian master equation (158!) (left side) and for (!59ll (right side). 




Figure 8. On the left, the fidehty between the simulated state pg and that according 
to the Markovian master equation ((58)) . The analog using the Markovian master 
equation (j59l) is plotted on the right. In both plots the parameters and legends are the 
same. 



From these results one concludes that when modeling a system with multiple 
baths at different temperatures equations (155]) and (151?]) are each accurate in their 
theoretically applicable regimes. However, for baths at the same temperature, it seems 
both equations give good results. A natural, and important, question is to ask is whether 
an intermediate range of couplings exist, such that neither ( 158]) or (159]) give useful results. 
In figure [9] the fidelity between the simulation and the Markovian master equation 
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states have been plotted for both equations at fixed time t = 100 as a function of the 
intercouphng strength (3. 
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Figure 9. Fidelity between the simulated state p'g'' and p^™'' according to the 
Markovian master equations (|58l) and (|59p at fixed time as a function of the coupling 
between the damped oscillators. 

We see that for the parameters shown on the plot, there is a small range between 
P ~ 0.01 — 0.02 where neither Markovian master equation obtains a high precision. 
However, note that this range becomes smaller as the coupling with the bath decreases, 
and so generally both master equations cover a good range of values of /3. 

3.3.1. Baths with the same temperature We now examine the role of the bath 
temperatures in more detail. Since the simulations seem to produce good results for 
both Markovian master equations when the temperature of the local baths are the same, 
regardless of the strength of the intercouphng, it is worth looking at why this happens. 
In the case of equation ( l59i) it is reasonable to expect that this will remain valid for 
small /9, because when /3 this equation approaches (!58l) if the bath temperatures and 
spectral densities are the same. That is, the off-diagonal terms of the matrices K^^'' and 
K^'^'^ do not contribute much, /3 ~ /3 and the rest of coefficients become approximately 
equal to those in fl58l ) Note this only happens under these conditions. 

Essentially the same argument applies to equation fISS]) in the large /3 limit. On the 
one hand, for a relatively small value of /3 (= 0.1) in comparison to w, the off-diagonal 
elements of the matrices K^^'^ and K'^^'^ in the master equation (1591) are unimportant in 
comparison with the diagonals. On the other hand, the diagonal terms are also alike for 
the same reason, and so both master equations will be quite similar. However note that 
at later times the behaviour of both equations start to differ, and the steady states are 
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not the same. By construction, the steady state of equation (1^ is the thermal state 
of the composed system [3l H] , whereas that of master equation (l58l) is not (ahhough it 
tends to the thermal state as /3 — )■ of course). Surprisingly the divergences between 
both equations, even for large times, are actually very small, see figure [TOl In some 
cases, while the steady state of (l58ll is not strictly thermal, the fidelity with that of (l59l) 
is more than 99.999%. 




Parameters: 



- Tsi = Ts2 = 20. 

- Tbi = Tb2 = 10. 
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Figure 10. Fidelity between states p^"^'' and Pg™^'' corresponding to Markovian 
master equations ([55|) and respectively. 

4. Driven Damped Harmonic Oscillator 

One situation which is also interesting to analyze is that of adding a driving term in 
the Hamiltonian of the damped oscillator. At this stage we consider again one single 
oscillator, damped by a thermal bath and driven by a coherent field (figure [TT]) . This is 
described by a semiclassical Hamiltonian in the rotating wave approximation: 

M M 

H{t) = Qa^a + ria^e''"^^* + ae'"^^^) + ^ Wj-ajaj + '^gj{a^aj + aa]), (64) 
here ul is the frequency of the incident field and r the Rabi frequency. 
Exact Solution 

To obtain the exact solution of this system let us consider for a moment the Schrodinger 
picture. 
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Figure 11. A single damped oscillator interacting with a classical incident field with 
Rabi frequency r. 



We solve this equation by means of the unitary transformation f/rot(i) = e where 

\ Making the substitution Itpit)) = f/rot(i)|^(i)) we 



immediately obtain 

d\m) 



dt 



i[Hrot - u,ot{t)H{t)uUt)]m)) = -^Ho\m), 



where Hq = {Q - ujL)a)a + r(a + a^) + Y^f^^ioOj - UL)a]aj + J2jLi dji^^^^j + ^^]) is 
time- independent. Returning to the Schrodinger picture, the evolution of the states is 
then, 

In order to avoid differential equations with time- dependent coefficients, we can study 
the evolution in a X-P time rotating frame; in that frame the annihilation (and creation) 
operators d = e~^^'°^^ae^^^°^^ will evolve according to 



a{t) = U^t,0)t 



Uit,0) 



^iHot^^-iHot_ 



That is 



M 

id = [d, Hq] = {Q — ujLjd + Qjdj + 
idj = [dj, Hq] = {uj — (jJLjdj + Qjd, 



(65) 
(66) 

which is quite similar to fl5^ but with the additional time- independent term r. Following 
the notation of section 12.11 we can write 

iA = WqA + 6, 

here b = (r, 0, . . . , 0)""" and Wq is found from (135!) as W — wlI. The solution of this 
system of differential equations is 



A{t) 



-iWot 



A{0) - i I dse 





iWosi 



Markovian Master Equations: A Critical Study 26 

If Wq is invertible this equation can be written as 

A{t) = e"^^"* [A{0) + Wo^b] - W^^b, (67) 

Analogously to ( l38l) and ( l39i) we find 

X{t) = T^X{0) - T^P{0) + T^Wo^b - (68) 
P{t) = T^X{0) + r]|P(0) + T^Wo\ (69) 

where and T° are as in f l40|) for W^o- Thus, by writing 

[t? n J' [ T^W.'b 

we find that the position and momentum expectation values evolve as 

R{t) = M°R{0) + B. (70) 

Note that in this case the first moments of the state change, despite (-R(O)) = 0. To 
calculate the evolution of the covariance matrix, we proceed in the same way as before, 

+B,J2KMm + l3.13„ (71) 

k 

and analogously for the solutions for {Rj{t)Ri{t)) and {Ri(t)){Rj(t)). Combining these 
terms, we find the B cancel and so, in a similar fashion to (144|) . (145!) and (146!) . 

Ci,i(t) = (A^?,C(0)Al?), 

Cl,M+2(t) = CM+2At) = iM'„CM%^,) 

CM+2,M+2it) = CMl,+,), (72) 

where, of course, A^? and as in fllSl) for Ai*^. 

4-2. Markovian Master Equations 

In order to derive a Markovian master equation for this system we must take account 
of two important details. First, since the Hamiltonian is time-dependent the generator 
of the master equation must also be time- dependent, 

whose solution defines a family of propagators S{t2,ti) such that 
Psih) = £it2,ti)psiti), 

£{t:,M)=£{hM)£{t2M)- 

These can be written formally as a time-ordered series 
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where T is the well-known time-ordering operator. Similarly to the case of time- 
independent equations it can be shown that the family S{t2,ti) is completely positive 
for all (^2 > ^i) if and only if £t has the Kossakowski-Lindblad form for any time t [T3] . 

The second problem is that there is an absence of rigorous methods to arrive at 
a completely positive master equation in the Markovian limit when the Hamiltonian 
is time-dependent, with the exception of adiabatic regimes of external perturbations 
[37] . Fortunately in this case, due to the simple periodic time-dependence of the 
Hamiltonian, we will be able to obtain Markovian master equations valid for large (to 
some degree) Rabi frequencies, even though the complexity of the problem has increased. 
In our derivation, we will distinguish between three cases: these will be when the Rabi 
frequency is very small; when the driving is far off resonance (la;^, — fi| ^0) and finally 
the identical case without the secular approximation. 

The details of the derivation are left for the Appendix B.3 but in these three cases 
we find a Markovian master equation with the structure 
d 



dt 



ps = -i[ila'^a + fe"^^^a + f*e 



where V is given by f l57|) . = + A is the same as for a single damped oscillator, and f 
is a renormalized Rabi frequency due to the effect of the bath. Note that as the incident 
field alters the position operator of the oscillator, which in turn couples to the bath, one 
should expect that the field is itself also effected by the environment. For small Rabi 
frequencies an argument similar to section 13.2.11 gives simply 

f = r, (73) 

whereas, when the driving field is far from resonance, \ujl — f2| ^ 0, we obtain 

A(fi) + i7((])' 



^+ o 

\l — ujl 

Finally, if we neglect the secular approximation, this regime yields 



(74) 



(75) 



Without entering into the details of the derivation, one sees that equations fl74|) and 



( !75l) are problematic on resonance — w^^l ~ 0. This is due to two approximations, one 
is the secular approximation in ( 1741) . and the other is the second order in the perturbative 
series. In the derivation in Appendix B.3 it is clear why in this case the series diverges 
for \VL — uA {). 



4-3. Study of the Approximations 

Note that in this case the range of validity of each equation is now more ambiguous than 
in previous sections where we have dealt with undriven systems. Which one is more 
appropriate is going to be discovered by simulation, although one could suppose that 
the more elaborate equations (ITill and ( |75l) would provide the better approximation. 



Markovian Master Equations: A Critical Study 



28 



However, there is still the question of how effective they are, and whether the additional 
effort required to obtain them is worthwhile in comparison to the simpler equation (I73|) . 

In addition note that in every case the covariance matrix is unaffected by the driving 
term, which only produce a change in the first moments. Furthermore, as the fidelity is 
invariant under unitary operations, we are always free to work in the frame rotating with 
the field. Therefore, all calculations can be performed with the rotating observables. 




Time Time 



Figure 12. Fidelity between pg and pg for different renormalized Rabi frequencies 
([73)l . (|74l) and ((75|) . An example of off resonance is shown on the left, whereas the plot 
on the right is close to resonance. 

In figure [12] the fidelities are plotted for close to and far from resonance. Compare 
the amount of disagreement with the fidelity of a single damped oscillator in figure El For 
global features, the more elaborate equation f l75|) works better in both cases, although 
the difference with fl73l) is very small. As expected, the choice of fTM]) is preferable 
to the choice of f l75]) when out of resonance, but gives quite poor results when close 
to resonance. However, when off resonance the difference among the three choices is 
essentially small. 

Given these results, it is worthwhile to look at how the fidelities at one fixed time 
vary as a function of the detunning, this is done in figure [13] (note we choose a large 
value for the time, so we avoid the potentially confusing effect due to the oscillatory 
behaviour depicted in figure [T2|) . 

Here we see that both (171]) and fail close to resonance, as was expected from 
the perturbative approach. Equation f[73|) gives good results due to the small Rabi 
frequency, however note in comparison to f l75l) the accuracy quickly drops off as we 
move away from ul — Q = 0. A similar effect can be seen when compared to ([74]) for 
larger detunnings. 

Finally, in figure [14] we test the dependency of the fidelities on the strength of the 
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Figure 13. Fidelity between pg and pg for different renormalized Rabi frequencies 
(|75)) . (17^ and ([75)1 as a function of tlic detunning. 



Rabi frequencies far from resonance. Here the worst behaviour is observed for f l73|) . as 
expected. 




(s) (m) 

Figure 14. Fidelity between pg and pg for different renormalized Rabi frequencies 
([73)) . d74l) and (|75)) as a function of the Rabi frequency. 



In summary, for the case of a driven damped harmonic oscillator the difference in 
accuracy among Markovian master equations is generally small. Equations f l7^ and fl75|) 
work better except in the case of resonance, where (I73|) gives more accurate results, as 
long as the Rabi frequency is small. The justification to use one equation over another 
will depend on the context and the accuracy which one wants to obtain, but given that 
the differences are so small the simplest choice fl73p seems to be the more "economical" 
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way to describe the dynamics. 
5. Conclusions 

We have obtained and studied the range of vahdity of different Markovian master 
equations for harmonic oscillators by means of exactly simulating the dynamics, and 
comparing the predictions with those obtained from evolving the system using the 
master equations. In particular, 

• We have clarified the possible detrimental effect of low temperatures on the 
Markovian treatment of a damped oscillator, showing that the Markovian master 
equation provides good accuracy regardless of the temperature of the bath. 

• We have shown that the system-environment state factorization assumption for all 
times has to be understood in general as an effective model by deriving the same 
equation using the projection operator technique. 

• We analysed two strategies for finding completely positive Markovian master 
equations for two harmonic oscillators coupled together under the effect of local 
baths, indicating that both are complementary in their range of validity. Moreover, 
when the temperature of the local baths is the same the difference between them 
is quite small. 

• In the same spirit, we derived time inhomogeneous completely positive Markovian 
master equations for a damped oscillator which is driven by an external semi- 
classical field. We studied the validity of each one and pointed out that completely 
positive dynamics can be obtained even without secular approximation (for these 
kinds of inhomogeneous equations). 

Despite the fact that we have focused on harmonic oscillator systems, the proposed 
method is general and we expect that non-harmonic systems should behave in a similar 
manner with respect to the validity of the equations. This suggest that the general 
conclusions made here are widely applicable to any other settings involving a weak 
interaction with an environment. 

In this regard, we hope that the present study may help in providing a better 
understanding and a transparent description of noise in interacting systems, including 
those situations where the strength of the internal system interaction is large. There 
are currently many quantum scenarios open to the use of these techniques, including 
realizations of harmonic and spin chains in systems of trapped ions [38] , superconducting 
qubits [39] and nitrogen- vacancy (NV) defects in diamond [ID]. 

Moreover, interacting systems subject to local reservoirs have been recently 
treated under the assumption of weak internal system interaction in theoretical studies 
ranging from the excitation transport properties of biomolecules [31] to the stability of 
topological codes for quantum information [32] • 
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Appendix A. Details of the simulation 

In order to make an appropriate comparison between the exact evolutions, such as those 
in sections 12. ![ 13.11 and 14. H and the corresponding master equations, we must make a 
careful choice of a number numerical parameters. In practice, however, this is not a 
difficult issue. The essential ingredient is to choose the couplings to the bath according 
to the desired spectral density. Throughout this paper, we have made the choice (13T|) . 

j 

The first step in picking gj is to remove the Dirac delta functions by integrating over a 
frequency range bounded by a frequency cut-off Wmax, 

J 

which means 

g], ^ aufe-'^^'/^'^^Auf 

due to the decomposition of the integral in terms of Riemann sums. We should also 
take care to set the range of oscillators, Wmax, large enough to cover fj3T|) significantly. 
For example, if we take ui = c, with c small, then one possible convention is to take 
Wmax such that J(wmax) = J{c), and so we neglect all possible oscillators with coupling 
constant less than a/ J(c) Awi. Another polisher convention is to take ooi and Wmax such 
that 

/ J{u)du = / J{u)du 

Jo J O^niax 

However, in practice this choice is not really a crucial point. 



Appendix B. Derivation of Markovian Master equations 

Appendix B.l. Two coupled damped harmonic oscillators, small fi 



We can derive Markovian master equations like fl58|) from the microscopic model by the 
following procedure. The Von Neumann equation in the interaction picture with respect 
to the free Hamiltonian iffree = -f^oi + + Hbi + Hb2 is 

j^p{t) = -^P[Vr2{t),pit)]-^o^\VsB{t),pit)] 

= /3Vi2(t)/^(t) + «V5B(t)p(t), (B.l) 
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where Yssit) = ViBi{t) + V2B2(t) and for simplicity we have assumed that the strength 
of the couphng to each bath is identical (the reader will note afterwards that this is not 
a crucial assumption). We now define the projector Vp(t) = Tt(^bi,B2) [p(t)] ® Pthi ® Pth2, 
along with Q = 1 —V. The application of the projection operators on ( IB.ip yields 

^Vp{t) = pVVi2{t)pit) + aVVsB{t)p{t), (B.2) 

j^Qm = /3QVi2(t)p(t) + aQVsB{t)m, (B.3) 
and so (c.f. section ll.ip we find a formal solution to the second equation as 

Qp{t) = g{t,to)Qp{to) + 13 [ dsg{t,s)QVu{s)Vp{s) 

J to 

+a [ dsg{t,s)QVsB{s)Vp{s), (B.4) 

J to 

where 

g(t,s) = 7-e^s'^*'s[/3Vi2(t')+"VsB(t')], 

Now the procedure is as follows, we introduce the identity 1 = P + Q in the second 
term of equation ( IB. 21) , 

j.Vp{t) = PVVu{t)m + aVVsB{t)Vp{t) + aVVsB{t)Qp{t), 

and insert the formal solution (IB. 40 into the last term. Recalling the condition ([2]) 
VVV = and again assuming an initial factorized state (Qp(to) = 0) we find 

^Vp(t) = /3VVu{t)m+ [ dslCi{t,s)Vp{s)+ [ dslC2{t,s)Vp{s), 

Jto Jto 

where here the kernels are 

IC^{t,s) = aPVVsB{t)g{t, s)QVu{s)V = 0, 
IC2{t,s) = a^VVsB{t)g{t,s)QVsB{s)V. 

The first vanishing because Vi2(s) commutes with V and QV = 0. If we consider the 
second kernel, weak coupling implies a > and so to second order in a and /3 this 
becomes 

/C2(t, s) = a'VVsB{t)QVsB{s)V + a'/S), 
which has exactly the same form as ([6]) and therefore the equation of motion becomes 

jVpit) = l3VVi2{t)p{t) + a" dsVVsB{t)VsB{s)Vp{s). 

Finally we note that 



Tr 



B1,B2 



ViBl(t)V2B2(t') (pthl (S) Pth2) 

TrBi[VlBi(t)pthi] TrB2[V'2B2(t')Pth2] = 0, 
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because our interactions individually hold TrBi[ViBipthi] = '~^T^B2[V2B2Pth2] = 0, so 
VV1B1V2B2V = VV2B2V1B1V = and then 

^Vp{t) = pVVu{t)m + ^ dsVViBi{t)ViBi{s)Vp{s) 



W [ dsVV2B2{t)V2B2is)Vp{s), 
J to 

which may be rewritten as 



j^psit) = -l[Vi2{t),Ps{t)] - I dt'TT^i[ViBl{t), [ViBlit'), [P5(t') ® Pthl(t')]] 



to 



- [ dt'TTB2[V2B2{t), [V2B2{t'), [ps{t') ® PtUt')]]- (B.5) 
J to 

The last quantity in the above equation is just a sum of the individual terms for 
each bath, which lead, under the standard procedure of section ll.H to the (interaction 
picture) local dissipators Vi and V2 and shifts of (l58l) . 

Appendix B.2. Two coupled damped harmonic oscillators, large (3 

First, let us write the Hamiltonian of the two oscillator system in a more convenient 
way 

We can diagonalize this quadratic form by means of a rotation to get 

H12 = n+b\h + n^blb2, 

where 



{n, + ^2) ± v/4/32 + {Q, - Q,)^ 
o±- , 

and the creation and annihilation operators in the rotated frame are given by 

hi = ai cos(q;) — 02 sin(Q;), 
62 = cti sin(a) + 02 cos(q;), 

with the angle specified by 
tan(aj 



{Q, - Q2) - ^/WTW^^W' 

The new operators satisfy the standard bosonic commutation rules [bi, bj] = 6ij, and so 
this is nothing more than the decomposition of an oscillatory system in normal modes. 
For simplicity, let us now take Qi = Q2 = ^, and so 

{bi = -^(cLi + 02) 
b2 = -02) 



Markovian Master Equations: A Critical Study 34 

note that RWA approximation implies Q ^ P so both normal mode frequencies are 
positive. 

We can reexpress the interactions with the baths in terms of these new operators, 

M 



^1^1 = E 7f [(^1 + + (^1 + &2)al,,, 

j = l V 

M 

92j 



^2B2 = E 71 [(^1 - ^2)«2, + ih - b,)ay, 

j = l V 

the benefit of this is that it allows us to easily deal with the interaction picture with 
respect to Hq = H12 + Hbi + Hb2- By following the method of section [LT] we obtain 
the analog of (|H]), 

d ^ ^ 

-Psit) = - J dt'TTBi[ViBi{t),[ViBi{t-s),ps{t)^ptu]] 

POO 

- dsTlB2[V2B2it),[V2B2it-s),psit)^ptt2]], (B.6) 

^0 

where we have noted VVibiV2B2'P = 'PV2B2^ibi'P = 0. Each of the above terms 
correspond, essentially, to one of a pair of two free harmonic oscillators with frequencies 
and , coupled to a common bath. Consequently, we can deal with them separately. 
Starting with the first term 

Ci{ps) = -j dsi:iBi[ViBi{t),[ViBi{t-s),ps{t-s)®pt^i]l (B.7) 
Jo 

we decompose the interaction in to eigenoperators of [-^12, ■] (see ( fTTi) ) 

^iBi = 5^A®5,, (B.8) 



with 



Al = 4=(&l+&2), A2 = ^{h\+hl), 



M M 

Bi = ^gija\,j, B2 = ^gijaij. (B.9) 
i=i i=i 

Notice the Ai operator can be written as Ai = Ai{il^)+Ai{il_), where Ai{il^) = bi/\/2 
and = 62/ are already the eigenoperators of [-^12, ■] with eigenvalues — i7+ 

and — respectively. Similarly A2 = A2{—fl+) + A2{—^-), with A2{—fl+) = b\/V^ 
and ^2(-^^-) = 4/^2, and so we can write flB.81) as 

ViBi = Y,Ak®Bu = Y, Ak{v) ®Bk = Y^ Al{v) ® bI (B.IO) 

k v,k u,k 

which in interaction picture becomes 

ViBiit) = Y^f^-'^'M^) ® Bk{t) = 5^e*'^*4(z/) ® Blit). 
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Now, for the first element of (1131) we have 

= Y.3l dse'^''-^^^>[nr{oo,,) + ll (B.ll) 

where the mean number of quanta in the first bath ni{uij) with frequency uij, is given 
by the Bose- Einstein distribution (151 p . Going to the continuous hmit we take M — )■ cxd 
and introduce the spectral density of the first bath Ji(a;) = Yljdij^i'^ ~ ^ij)^ 



ri,i(z/)=/ dcoJiiu) rfse'(""")^[ni(u;) + l]. 
Jo Jo 

Now using the well-know formula from distribution theory, 

^ dxe'^'y = 7i6{y) + iP.V. , 

and assuming i/ > 0, we split into real and imaginary parts, 
ri,i(z/) = liiiy)[ni{'') + 1] + ^[Ai(z/) + A;(z/)], 



where 



7i(z/) = 7rJi(z/), 



u — u 

00 



a;(.) = p.v. / dc^^^MM^). (B.12) 



U — 00 

Similar calculations give {u > 0) 

ri,2(-/^) = r2,i(z/) = 0, (B.13) 

r2,2(-/^) = 7i('^)^i(^^) - ^a;(z/). (B.14) 

Thus, equation flB.7p becomes 

+e^('^'-^)*r2,2(/^)[A2(z/)p5(t),4(^^')] 

+e^('^-'^')*q2(^)[^2(^'),P5(t)4(^)]- (B.15) 

Next we perform the secular approximation; the cross terms u' ^ u in the above 

expression, which go as e^^''**, can be neglected provided that 2/3 is large in comparison 
with the inverse of the relaxation rate {(3 ^ a) and so we obtain 

A(P5) = -^^^[b\b^, psit)] - ^^l^[blh, Mt)] 

+7i(fi+)[ni(fi+) + 1] Upsit)b\ - l{b%,psit)} 
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+7i(fi+)ni(fi+) {h\ps{t)h - \{hihlps{t)} 
+7l(^]_)[nl(fi_) + 1] (h2Ps{t)h\ - \{blh2rpsit)} 



2 



+7i(l^-)ni(fi_) Ups{t)h2 - \{h2hlps{t)}\ . 



Returning to equation flB.6p . for the second term, 



/■oo 

i^2{ps) = - dsTTBl[V2B2{t), [V2B2{t - s) , ps{t) ®Pth2]], 

Jo 

the situation is essentially the same, since the minus sign in 62 only modifies the cross 
terms, which we neglect in the secular approximation. Following similar steps as in the 
above we obtain the same form fIB.lGP for C2, with the replacements 71 — 72, Ai — t- A2 
and fii — )■ n2, where the subscript 2 refers to the corresponding expression with the 
spectral density and temperature of the second bath. Therefore putting together both 
quantities, and returning to the Schrodinger picture 

j^psit) = -I [n^ + Ai(l]+)/2 + A2(l]+)/2] [h\h,ps{t)] 

-I [^2 + Ai(f]_)/2 + A2(f]-)/2] [6^62, ps{t)] 

+{7l(^]+)[r^l(^]+) + 1] +72(fi+)[n2(^^+) + 1]} Ups{t)h\ - ]-{h\h,, psit)} 



+ [7i(l^+)ni(fi+) +72(^^+)r22(^^+)] {h\ps{t)hi - \{hhlPsit)} 
+{7l(^_)[r^l(^]_) + 1] +72(^]-)[n2(^^_) + 1]} {h2Ps{t)h\ - \{h\h2,ps{t)} 

+ [7i(fi_)ni(fi_) +72(r]_)n2(^^-)] {h\ps{t)h2 - \{h2hlps{t)}^ ■ (B.16) 

It is manifestly clear that this equation is of the Kossakowski-Lindblad form. Finally, 
we rewrite the operators hi and 62 in terms of ai and 02 to arrive at equation (!59|) . 

It is worth mentioning that similar equations for coupled harmonic oscillators have 
been given previously (see for example |l3l HI]), but not in the Kossakowski-Lindblad 
form, since in those derivations the secular approximation is not taken. 

Appendix B.3. Driven damped harmonic oscillator 

To derive a completely positive Markovian master equation valid for large Rabi 
frequencies r we must work in the interaction picture generated by the unitary 
propagator U{ti,to) = Te^^^^o ^^^^ , where 

M 

Hiit) = Qa^a + r{a^ e'^^^ + ae''^^*) + ^w^aja^. (B.17) 
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Taking to = without lost of generality, the time-evolution equation for p(t) = 
W{t,0)p{t)U{t,0) is 

'pit) = -^[V{t)rp{t)l (B.18) 

so by following the analogous procedure for time- independent generators, one 
immediately deals with the problem that is not clear whether there exists a similar 
eigenoperator decomposition for V{t) = W (t,0)VU{t,0) {V = JZjLidji^^'^C'j + 0.(1])) as 
in ( ITOl) and ( ITT]) . Note however that the operator Ai{t) = d{t) satisfies a differential 
equation with periodic terms 

ia{t) = [a{t), Ho{t)] = m{t) + re-'^^^K (B.19) 

This kind of equation can be studied with the well-established Floquet theory (see for 
example [261 |2l]), particularly it is possible to predict if its solution is a periodic function. 
In such a case, the operator in the new picture would have a formal decomposition 
similar to that in ( ITOj) and ( ITTi) . such that Akit) = Ylu '^kiy)^^^^^ ■ This would then 
allow us to follow a similar procedure to that for time-independent Hamiltonians. Note 
that the importance of such a decomposition is that the operators Ak^u) are themselves 
time- independent. Such ideas have already been used before in, for instance, ^5} B6]. 

The solution to equation ( IB. 191) . with the initial condition 5(0) = a and for Vt ^ ul 
is given by 

a{t) = -—^ , (B.20) 

so in this case the solution is periodic and the desired decomposition Ai(t) = 

T.^M^Y"' is 

A^{t) = A^{ujL)e'''^'' + A2{n)e-'''\ 

where Ai{uji) = and Ai{Q) = a — ^^'^^ 1 = a — Ai{ujl). Similarly 

A2{t) = A2i-UL)e'^^' + A2(-fi)e^™, 

with A2(— wl) = j^T^l = Ai{ijJl) and A2{—Q) = a'' — j^^l = a — A2{—col)- Thus we 
get an equation analogous to ( IB.lSp . where the coefficients are: 

Tniiy) = l{iy)[n{iy) + l] + i[A{u) + A'{u)], {v > 0) 
Tu{i^) = T2i{u) = 0, 

T22ii^) =l{-i^)ni~u) -iA\-i^) (z/<0). 

Before continuing note that in the perturbative series of ( IB. 181) . the "strength" of the 
interaction V{t) is now not solely dependent on the coupling with the bath. This is 
because the operators A^u) depend linearly on , so when this ratio becomes large 
we expect that the approximation breaks down, i.e. for r 3> 1 or very close to resonance 
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Next we assume that the detunning is large enough Icj^ — f2| ^ a, \ul — Q\'^ ^ ar 
in order to make the secular approximation and after some tedious, but straightforward, 
algebra we find the master equation in the interaction picture to be 



d ^ 
dtp' 



-i[A(fi)afa 



a + a 



7(f2)r a — a"!" 



Ps] + D{~ps 



(B.21) 



(jJl — Vt i 

where T>{-) has again the form of (J571) . Finally, on returning to the Schrodinger picture 
we have, 

'^ps = -z[H,{t),ps] + U{t,OypsU\t,0) 



dt 



Dips 



where f2 = i7 + A{Q) and 



1 + 



A{n) + ij{n) 



(B.22) 



(B.23) 



So in this master equation the Rabi frequency is renormalized by the effect of the bath. 
It is worth noting that at first order in r and the coupling a we obtain equation f l73|) . 
This is as expected, given the arguments in section 13.2. II 

For an arbitrary driving frequency a Markovian master equation is difficult to obtain 
as we cannot, in general, make the secular approximation (apart from the perturbative 
condition — Q\ ^ 0). This can be illustrated in the extreme case of resonance 
ujl = Solving equation ( 1B.19P under this condition we find 

d{t)=e-'^\a-irt), (B.24) 

and so one can see that d{t) is not a periodic function, so the desired decomposition 
as a sum of exponentials with time- independent coefficients does not exist. On the 
other hand, the decomposition (IB.20p tends to flB.24p in the limit ul — )■ f2, so we may 
attempt to work with this decomposition and wonder whether on resonance the new 
master equation holds in this limit as well (in fact, we have shown that this is not true 
in section . The only problem to deal with is the possible lack of positivity due to 
the absence of the secular approximation. However, note that in this particular case 
only a commutator term arises from the cross terms in the analog of equation fIB.lSp . 
so positivity is not lost. In fact, we obtain an equation similar to ( IB. 221) except for an 
additional correction to the Rabi frequency: 

AiSl) + i-iiSl) A(wi)+i7(wi) 



Note that to first order in r and a we again obtain the equation ( |56 



(B.25) 
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